function h = lpcvd_mode2h_2D_sincos(z1,z2,z3,z4,L)
h = zeros(L,L);
MaxMode = size(z1,1)-1;
dx = pi/L;
x = 0:dx:pi-dx;
y = 0:dx:pi-dx;
for i = 1:L
    for j = 1:L
        h(i,j) = 0.0;
        for m = 0:MaxMode
            for n = 0:MaxMode
                h(i,j) = h(i,j)+z1(m+1,n+1)*fun_phi1mn(m,n,x(i),y(j),L)...
                         +z2(m+1,n+1)*fun_phi2mn(m,n,x(i),y(j),L)...
                         +z3(m+1,n+1)*fun_phi3mn(m,n,x(i),y(j),L)...
                         +z4(m+1,n+1)*fun_phi4mn(m,n,x(i),y(j),L);
            end
        end
    end
end
%------------------------------------
    function phi1 = fun_phi1mn(m,n,x,y,L)
        phi1 = (2/L)*sin(2*m*x)*(sin(2*n*y))';
    end

    function phi2 = fun_phi2mn(m,n,x,y,L)
        if (m== 0 || m == L/2) && (n==0||n==L/2)
            coeff = 1/L;
        elseif m~=0 && m~=L/2 && n~=0 && n~=L/2
            coeff = 2/L;            
        else
            coeff = sqrt(2)/L;
        end
        phi2 = coeff*cos(2*m*x)*(cos(2*n*y))';
    end

    function phi3 = fun_phi3mn(m,n,x,y,L)
        if n == 0 || n == L/2
            coeff = sqrt(2)/L;
        else
            coeff = 2/L;
        end
        phi3 = coeff*sin(2*m*x)*(cos(2*n*y))';
    end

    function phi4 = fun_phi4mn(m,n,x,y,L)
        if m ==0 || m == L/2
            coeff = sqrt(2)/L;
        else
            coeff = 2/L;
        end
        phi4 = coeff*cos(2*m*x)*(sin(2*n*y))';
    end
end